clear;clc;
tic;

%% 载入图案数据
img = importdata('img.mat'); % 载入图案

for n = 1:2
    %% 调用CSTStudio并初始化
    A = img(:, :, n); % 下层图案
    a = img(:, :, n + 2); % 上层图案
    path = pwd; % 获取当前m文件夹路径
    cst = actxserver('CSTStudio.application');
    mws = invoke(cst, 'NewMWS');
    app = invoke(mws, 'GetApplicationName');
    ver = invoke(mws, 'GetApplicationVersion');
    invoke(mws, 'FileNew');
    invoke(mws, 'DeleteResults'); % 删除之前的结果，防止程序停止运行

    %% 仿真参数设置
    Px = 64; %结构X向边长
    Py = 64; %结构Y向边长
    h = 2; %基地层金属厚度
    h1 = 0.5; % 顶层金属层厚度
    h2 = 0.5; %底层金属层厚度
    l1 = 35; %顶层金属边长
    l2 = 50; %底层金属边长
    w = 5; %金属宽度
    t1 = 5; %顶层介质层厚度
    t2 = 5; %底层介质层厚度
    Frq = [2.5, 3.5]; % 工作频率
    little_length = 1; % 图案小块边长

    % 在CST中加入结构参数，方便后续查看
    invoke(mws, 'StoreParameter', 'Px', Px);
    invoke(mws, 'StoreParameter', 'Py', Py);
    invoke(mws, 'StoreParameter', 'h', h);
    invoke(mws, 'StoreParameter', 'h1', h1);
    invoke(mws, 'StoreParameter', 'h2', h2);
    invoke(mws, 'StoreParameter', 'l1', l1);
    invoke(mws, 'StoreParameter', 'l2', l2);
    invoke(mws, 'StoreParameter', 'w', w);
    invoke(mws, 'StoreParameter', 't1', t1);
    invoke(mws, 'StoreParameter', 't2', t2);

    %% 单位格式设置
    sCommand = '';
    sCommand = [sCommand 'With Units'];
    sCommand = [sCommand 10 '.Geometry "um"'];
    sCommand = [sCommand 10 '.Frequency "THz"'];
    sCommand = [sCommand 10 '.Time "ns"'];
    sCommand = [sCommand 10 '.TemperatureUnit "kelvin"'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define Units', sCommand);

    %% 频率范围设置
    sCommand = '';
    sCommand = [sCommand 'Solver.FrequencyRange "',num2str(Frq(1)),'", "', num2str(Frq(2)), '"'];
    invoke(mws, 'AddToHistory', 'Define Frequency Range', sCommand);

    %% 材料设置
    % 1. 背景材料设置
    sCommand = '';
    sCommand = [sCommand 'With Background'];
    sCommand = [sCommand 10 '.ResetBackground'];
    sCommand = [sCommand 10 '.Type "Normal"'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define Background', sCommand);

    % 2. 材料1Gold设置
    sigma = 4.561e+007;
    sCommand = '';
    sCommand = [sCommand 'With Material'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Gold"'];
    sCommand = [sCommand 10 '.FrqType "all"'];
    sCommand = [sCommand 10 '.Type "Lossy metal"'];
    sCommand = [sCommand 10 '.Sigma "', num2str(sigma), '"'];
    sCommand = [sCommand 10 '.Mu "1.0"'];
    sCommand = [sCommand 10 '.Rho "19320"'];
    sCommand = [sCommand 10 '.ThermalConductivity "314"'];
    sCommand = [sCommand 10 '.HeatCapacity "0.13"'];
    sCommand = [sCommand 10 '.YoungsModulus "78"'];
    sCommand = [sCommand 10 '.PoissonsRatio "0.42"'];
    sCommand = [sCommand 10 '.ThermalExpansionRate "14"'];
    sCommand = [sCommand 10 '.Colour "1","1","0"'];
    sCommand = [sCommand 10 '.Create'];
    sCommand = [sCommand 10 '.Delete "Vacuum"'];
    sCommand = [sCommand 10 '.Delete "PEC"'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define Material:Gold', sCommand);

    % 3. 材料2Dielectric设置
    Mu = 2.0;
    sCommand = '';
    sCommand = [sCommand 'With Material'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Dielectric"'];
    sCommand = [sCommand 10 '.FrqType "all"'];
    sCommand = [sCommand 10 '.Type "Normal"'];
    sCommand = [sCommand 10 '.Epsilon "1.0"'];
    sCommand = [sCommand 10 '.Mu "', num2str(Mu), '"'];
    sCommand = [sCommand 10 '.Create'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define Material:Dielectric', sCommand);

    % 4. 划分组成
    sCommand = '';
    sCommand = [sCommand 'Component.New "Dielectric"'];
    invoke(mws, 'AddToHistory', 'New Component', sCommand);
    sCommand = '';
    sCommand = [sCommand 'Component.New "Plasmonic"'];
    invoke(mws, 'AddToHistory', 'New Component', sCommand);
    sCommand = '';
    sCommand = [sCommand 'Component.New "Substrate"'];
    invoke(mws, 'AddToHistory', 'New Component', sCommand);

    %% 建立模型
    %% 1. 建立底层金属模型
    sCommand = '';
    sCommand = [sCommand 'With Brick'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Substrate"'];
    sCommand = [sCommand 10 '.Component "Substrate"'];
    sCommand = [sCommand 10 '.Material "Gold"'];
    sCommand = [sCommand 10 '.Xrange ', '"0", "Px"'];
    sCommand = [sCommand 10 '.Yrange ', '"0", "Py"'];
    sCommand = [sCommand 10 '.Zrange ', '"-h", "0"'];
    sCommand = [sCommand 10 '.Create'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define brick: Substrate', sCommand);

    %% 2. 建立下层介质层模型
    sCommand = '';
    sCommand = [sCommand 'With Brick'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Bottom Dielectric"'];
    sCommand = [sCommand 10 '.Component "Dielectric"'];
    sCommand = [sCommand 10 '.Material "Dielectric"'];
    sCommand = [sCommand 10 '.Xrange ', '"0", "Px"'];
    sCommand = [sCommand 10 '.Yrange ', '"0", "Py"'];
    sCommand = [sCommand 10 '.Zrange ', '"0", "t2"'];
    sCommand = [sCommand 10 '.Create'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define brick: Bottom Dielectric', sCommand);

    %% 3.1 建立下层金属层
    k = 1; % 用于统计小块的个数

    for i = 1:64

        for j = 1:64

            if A(i, j) == 1
                sCommand = '';
                sCommand = [sCommand 'With Brick'];
                sCommand = [sCommand 10 '.Reset'];
                sCommand = [sCommand 10 '.Name ', ['"Ring', num2str(k), '"']];
                sCommand = [sCommand 10 '.Component "Plasmonic"'];
                sCommand = [sCommand 10 '.Material "Gold"'];
                sCommand = [sCommand 10 '.Xrange "',num2str(little_length*(i-1)),'", "', num2str(little_length * (i)), '"'];
                sCommand = [sCommand 10 '.Yrange "',num2str(little_length*(j-1)),'", "', num2str(little_length * (j)), '"'];
                sCommand = [sCommand 10 '.Zrange "t2","t2+h2"'];
                sCommand = [sCommand 10 '.Create'];
                sCommand = [sCommand 10 'End With'];
                invoke(mws, 'AddToHistory', ['Define Brick: Bottom Plasmonic', num2str(k)], sCommand);
                k = k + 1;
            end

        end

    end

    % 3.2 合并下层金属
    for i = 1:k - 2
        sCommand = '';
        sCommand = [sCommand 'Solid.Add "Plasmonic:Ring1", ', ['"Plasmonic:Ring', num2str(i + 1), '"']];
        invoke(mws, 'AddToHistory', 'Add Bottm Metallic', sCommand);
    end

    %% 4.1 创建上层介质层
    sCommand = '';
    sCommand = [sCommand 'With Brick'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Top Dielectric"'];
    sCommand = [sCommand 10 '.Component "Dielectric"'];
    sCommand = [sCommand 10 '.Material "Dielectric"'];
    sCommand = [sCommand 10 '.Xrange "0","Px"'];
    sCommand = [sCommand 10 '.Yrange "0","Py"'];
    sCommand = [sCommand 10 '.Zrange "t2","t2+t1"'];
    sCommand = [sCommand 10 '.Create'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Define brick: Top Dielectric', sCommand);

    % 3.2 在上层介质层中创建相同的下层金属层并除去
    sCommand = '';
    sCommand = [sCommand 'With Transform'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Name "Plasmonic:Ring1"'];
    sCommand = [sCommand 10 '.Vector "0","0","0"'];
    sCommand = [sCommand 10 '.UsePickedPoints "False"'];
    sCommand = [sCommand 10 '.InvertPickedPoints "False"'];
    sCommand = [sCommand 10 '.MultipleObjects "True"'];
    sCommand = [sCommand 10 '.GroupObjects "False"'];
    sCommand = [sCommand 10 '.Repetitions "1"'];
    sCommand = [sCommand 10 '.MultipleSelection "False"'];
    sCommand = [sCommand 10 '.Destination ""'];
    sCommand = [sCommand 10 '.Material ""'];
    sCommand = [sCommand 10 '.Transform "Shape", "Translate"'];
    sCommand = [sCommand 10 'End With']
    invoke(mws, 'AddToHistory', 'Paste Bottom Metallic', sCommand);

    sCommand = '';
    sCommand = [sCommand 'Solid.ChangeComponent "Plasmonic:Ring1", "Dielectric"'];
    invoke(mws, 'AddToHistory', 'Subtract Bottom Metallic', sCommand);
    sCommand = '';
    sCommand = [sCommand 'Solid.Subtract "Dielectric:Top Dielectric", "Dielectric:Ring1"'];
    invoke(mws, 'AddToHistory', 'Subtract Bottom Metallic', sCommand);

    %% 5.1 创建顶层金属层
    clear k;
    k = 1;

    for i = 1:64

        for j = 1:64

            if a(i, j) == 1
                sCommand = '';
                sCommand = [sCommand 'With Brick'];
                sCommand = [sCommand 10 '.Reset'];
                sCommand = [sCommand 10 '.Name ', ['"Top Metallic', num2str(k), '"']];
                sCommand = [sCommand 10 '.Component "Plasmonic"'];
                sCommand = [sCommand 10 '.Material "Gold"'];
                sCommand = [sCommand 10 '.Xrange "',num2str(little_length*(i-1)),'", "', num2str(little_length * (i)), '"'];
                sCommand = [sCommand 10 '.Yrange "',num2str(little_length*(j-1)),'", "', num2str(little_length * (j)), '"'];
                sCommand = [sCommand 10 '.Zrange "t1+t2","t1+t2+h1"'];
                sCommand = [sCommand 10 '.Create'];
                sCommand = [sCommand 10 'End With'];
                invoke(mws, 'AddToHistory', ['Define Brick: Top Plasmonic', num2str(k)], sCommand);
                k = k + 1;
            end

        end

    end

    % 5.2 合并顶层金属
    for i = 1:k - 2
        sCommand = '';
        sCommand = [sCommand 'Solid.Add "Plasmonic:Top Metallic1", ', ['"Plasmonic:Top Metallic', num2str(i + 1), '"']];
        invoke(mws, 'AddToHistory', 'Add Top Metallic', sCommand);
    end

    %% 定义边界条件及端口
    sCommand = '';
    sCommand = [sCommand 'With FloquetPort'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.Port "Zmin"'];
    sCommand = [sCommand 10 '.SetNumberOfModesConsidered "2"'];
    sCommand = [sCommand 10 '.SetUseCircularPolarization "True"'];
    sCommand = [sCommand 10 '.Port "Zmax"'];
    sCommand = [sCommand 10 '.SetNumberOfModesConsidered "2"'];
    sCommand = [sCommand 10 '.SetUseCircularPolarization "True"'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Set FloquetPort', sCommand);

    sCommand = '';
    sCommand = [sCommand 'With Boundary'];
    sCommand = [sCommand 10 '.Xmin "unit cell"'];
    sCommand = [sCommand 10 '.Xmax "unit cell"'];
    sCommand = [sCommand 10 '.Ymin "unit cell"'];
    sCommand = [sCommand 10 '.Ymax "unit cell"'];
    sCommand = [sCommand 10 '.Zmin "expanded open"'];
    sCommand = [sCommand 10 '.Zmax "expanded open"'];
    sCommand = [sCommand 10 'End With'];
    invoke(mws, 'AddToHistory', 'Set Boundary', sCommand);

    %% 设置求解器
    sCommand = '';
    sCommand = [sCommand 'With FDSolver'];
    sCommand = [sCommand 10 '.Reset'];
    sCommand = [sCommand 10 '.LimitIterations "True"'];
    sCommand = [sCommand 10 '.MaxIterations "10"'];
    sCommand = [sCommand 10 '.SetNumberOfResultDataSamples "1001"'];
    sCommand = [sCommand 10 '.SetResultDataSamplingMode "Frequency (linear)"'];
    sCommand = [sCommand 10 '.Start'];
    sCommand = [sCommand 10 'End With']
    invoke(mws, 'AddToHistory', 'Define FDSolver', sCommand)

    %% 后处理
    CST_Post = CST_MicrowaveStudio();
    [freq, sparam, stype] = CST_Post.getSParameters;
    % 得到的S参数依次为：
    % Szmax1,zmax1    Szmax2,zmax1    Szmin1,zmax1    Szmin2,max1
    % Szmax1,zmax2    Szmax2,zmax2    Szmin1,zmax2    Szmin2,zmax2
    % Szmax1,zmin1     Szmax2,zmin1     Szmin1,zmin1    Szmin2,zmin1
    % Szmax1,zmin2     Szmax2,zmin2    Szmin1,zmin2    Szmin2,zmin2
    % Szmax1,zmax1 == Szmax2,zmax2
    % 取CD = Szmax2,zmax1-Szmax1,zmax2
    % frequency = freq(:,1); % 获取频率点
    % save('frequency.mat','frequency');
    real = real(sparam); % 获取实部
    imag = imag(sparam); % 获取虚部
    linear = sqrt(real .* real + imag .* imag); % 计算linear值
    Szmax1zmax1 = linear(:, 1);
    Szmax2zmax1 = linear(:, 2);
    Szmax1zmax2 = linear(:, 5);
    CD = Szmax2zmax1 .* Szmax2zmax1 - Szmax1zmax2 .* Szmax1zmax2;
    real_name = [path, '\result\real\real_', num2str(n), '.mat'];
    imag_name = [path, '\result\imag\imag_', num2str(n), '.mat'];
    Szmax1zmax1_name = [path, '\result\S-Parameters\Szmax1zmax1_', num2str(n), '.mat'];
    Szmax2zmax1_name = [path, '\result\S-Parameters\Szmax2zmax1_', num2str(n), '.mat'];
    Szmax1zmax2_name = [path, '\result\S-Parameters\Szmax1zmax2_', num2str(n), '.mat'];
    CD_name = [path, '\result\CD\CD_', num2str(n), '.mat'];
    save(real_name, 'real');
    save(imag_name, 'imag');
    save(Szmax1zmax1_name, 'Szmax1zmax1');
    save(Szmax2zmax1_name, 'Szmax2zmax1');
    save(Szmax1zmax2_name, 'Szmax1zmax2');
    save(CD_name, 'CD')
    clear freq sparam stype real imag linear Szmax1zmax1 Szmax2zmax1 Szmax1zmax2 CD
    % 保存并退出
    % invoke(mws,'Save');
    filename = ['\metamaterials\metamaterial', num2str(n), '.cst']; % 新建的CST文件名字
    fullname = [path filename];
    invoke(mws, 'SaveAs', fullname, 'True'); % 保存到目前为止的结果
    invoke(mws, 'Quit');

end

toc
